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Abstract 


We  consider  the  steady  state  equations  for  a  compressible  fluid.  Since  we  wish  to  solve  for  a  range  of 
speeds  we  must  consider  the  equations  in  conservation  form.  For  transonic  speeds  these  equations  are  of 
mixed  type.  Hence,  the  usual  approach  is  to  add  time  derivatives  to  the  steady  state  equations  and  then 
march  these  equations  in  time.  One  then  adds  a  time  derivative  of  the  density  to  the  continuity  equation,  a 
derivative  of  the  momentum  to  the  momentum  equation  and  a  derivative  of  the  total  energy  to  the  energy 
equation.  This  choice  is  dictated  by  the  time  consistent  equations.  However,  since  we  are  only  interested 
in  the  steady  state  this  is  not  necessary.  Thus  we  shall  consider  the  possibilty  of  adding  a  time  derivative 
of  the  pressure  to  the  continuity  equation  and  similar  modifications  for  the  energy  equation.  This  can  then 
be  generalized  to  adding  combinations  of  time  derivatives  to  each  equation  since  these  vanish  in  the  steady 
state.  When  using  acceleration  techniques  such  as  residual  smoothing,  multigrid,  etc.  these  are  applied  to 
the  pressure  rather  than  the  density.  Hence,  the  code  duplicates  the  behavior  of  the  incompressible  equations 
for  low  speeds. 
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Introduction 


It  is  well  known,  that  it  is  difficult  to  solve  the  compressible  equations  for  low  Mach  numbers.  For  an 
explicit  scheme  this  is  easily  seen  by  inspecting  the  time  steps.  For  stability,  the  time  step  must  be  chosen 
inversely  proportional  to  the  largest  eigenvalue  of  the  system  which,  for  slow  flows,  is  approximately  the 
speed  of  sound,  c.  However,  other  waves  are  convected  at  the  fluid  speed,  u  ,  which  is  much  slower.  Hence, 
these  waves  don’t  change  very  much  over  a  time  step.  Thousands  of  time  steps  may  be  required  to  reach  a 
steady  state.  Should  one  try  a  multigrid  acceleration  one  finds  that  the  same  disparity  in  wave  speeds  slows 
down  the  multigrid  acceleration.  With  an  implicit  method  an  ADI  factorization  is  generally  used  so  that 
one  can  easily  invert  the  implicit  factors.  The  use  of  ADI  introduces  factorization  errors  which  again  slows 
down  the  convergence  rate  when  there  are  wave  speeds  of  very  different  magnitudes  [7]  . 

We  consider  systems  of  the  form 

+  /ar  +  5^?/  =  0. 

Our  analysis  will  be  based  on  the  linearized  equations  so  that  the  conservation  form  does  not  appear  in  the 
analysis  though  it  does  appear  in  the  final  numerical  approximation.  This  system  is  now  replaced  by 


P  +  /x  +ifj,  =  0, 


or  in  linearized  form 

+  Bwy  =  0,  (1) 

with  A  and  B  constant  matrices. 

In  order  for  this  system  to  be  equivalent  to  the  original  system,  in  the  steady  state,  we  demand  that 
have  an  inverse.  This  only  need  be  true  in  the  flow  regime  under  consideration.  We  shall  see  later  that 
frequently  P  is  singular  at  stagnation  points  and  also  along  sonic  lines.  Thus,  the  final  preconditioner  will 
be  smoothed  out  in  the  vicinity  of  points  where  M=0  or  M=l. 

Assuming  the  steady  state  has  a  unique  solution,  it  does  not  matter  which  system  we  march  to  a  steady 
state.  We  shall  later  see  that  for  the  finite  difference  approximations  the  steady  state  solutions  are  not 
necessarily  the  same  and  usually  the  preconditioned  system  leads  to  a  better  behaved  steady  state. 


The  Incompressible  Limit 

The  time  dependent  two  dimensional  Euler  equations  can  be  written  as 


Pi  +  upx  +  vpy  +  pa^{ux  +  =  0 

Px 

Ui  +  UUx  -f  VUy  +  —  =  0 

Py  ^ 

Vi  -f  UVx  +  VVy  +  ~  =  0 

Si  “h  uSx  vSy  —  0 


(2) 


The  form  of  this  system  is  unchanged  if  we  nondimensionalize  the  equations.  From  now  on  we  shall 
assume  that  u,v,p^p  are  nondimensional  quantities  where  the  dimensional  variables  are  nondimensionalized 
by  Following  [4]  we  define  e  =  ^  =  M*  .  If  the  fluid  is  isentropic  then 


(3) 


Hence,  as  e  goes  to  zero  the  speed  of  sound,  a,  goes  to  infinity  and  so  the  first  equation  in  (2)  reduces  to 

Ux  +  %  =  0. 
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It  was  pointed  out  in  ([10],  [11])  that  these  equations  can  be  symmetrized  by  using  ^  as  the  independent 
variable  rather  than  dp  .  Hence,  we  define  d>  hy  d<f)=  For  isentropic  flow  both  p  and  a  are  functions  only 

of  the  density  and  so  using  (3)  this  can  be  integrated  explicitly.  This  gives  <}>  =  .  As  the  Mach  number 

goes  to  zero  (f)  tends  to  infinity  and  therefore,  Gustafsson  and  Stoor  [4]  subtract  a  constant  and  define 


2  ^ 


(4) 


This  amounts  to  specifying  the  constant  in  the  integration  of  d<j)  from  dp.  They  then  prove,  using  energy 
methods,  that 

,  dpincompressible 

- di - 

Hence,  (j)  and  all  its  derivatives  behave  as  0{M)  as  M  0.  Since  p \  and  using  the  definition  of  d<f)  this 
is  equivalent  to 

dpcompressible  ^  dpin compressible  (5) 

We  now  consider  how  to  construct  a  matrix  artificial  viscosity  that  will  enable  us  to  reach  the  income 
pressible  limit.  Consider 


P  ^Wt  +  U  +  9y  =  h  [{QxWr)x  +  {Q2Wy)y] 


(6) 


We  wish  to  find  the  dependence  of  P  and  Qi  on  the  Mach  number  as  M  0  so  that  we  get  the  proper 
convergence.  We  therefore  consider  the  isentropic  equations  based  on  w  =  see  (4).  This  has  the 

symmetric  form 


P  ^Wt  + 


0>u 

U12 

diz  ^ 

^12 

^22 

0‘2S 

Wx 

an 

«23 

^33  y 

1 

+ 


bi2  hs 

bl2  ^22  ^23 

^11  hs  hs 


=  h[{QlW^)a:  -f  (<92^y)2/] 


As  M  0,ai2  and  613  =  0{l/M)  while  d(j)  =  0{M)  while  all  other  quantities  are  bounded.  Hence,  the 
leading  terms  in  the  first  equation  are  all  0(1/ M)  while  they  are  0(1)  for  the  second  and  third  equations. 
Multiplying  the  first  equation  by  M  and  taking  the  limit  we  get  u^:  +  Vy  for  the  space  derivatives  on  the  left 
hand  side.  Using  d<f)  =  0{M)y  du  =  0(1),  dt;  =  0(1)  we  see  that  a  necessary  condition  for  convergence  as 


(7) 


equation 

^’i  P fx  {\Q\'l^x^x 

Let  -4  =  1^.  Since  we  are  updating  Pf  we  should  have  Q  =  PA,  However,  this  is  not  in  conservation  form 
at  the  steady  state.  Instead  we  consider  artificial  viscosities  of  the  form 

P-^u,  +  U^{p-WPA\u^), 


M  — ^  0  is  that  P  have  the  form 

P-SQi,Q2~  0(4) 
V  0(4) 

0(4) 

Oiif) 

0(1) 

0(1) 

0(1) 

0(1) 

The  artificial  viscosity 

matrices 

Qi  are 

or 

Ui  +  PU  ~  P{P-\\PA\u^)^ 
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This  would  be  equivalent  to  the  original  form  if  P  were  constant.  Instead  we  have  terms  like  PjP  ^i+i/2 
that  appear. 

We  note  that  the  conditions  on  the  matrix  (P”^|PD|)comp  are  not  satisfied  by  the  non-preconditioned 
Roe  matrices.  Furthermore,  even  reasonable  preconditioners  need  not  satisfy  these  conditions.  Consider,  for 
example,  the  one  dimensional  system 

+  Awj;  =  h{Qwx)x 

A  reasonable  choice  is  P”^  =  \A\  i.e.  P  =  \A~^\.  In  this  case  all  the  wave  speeds  of  PA  are  ±1.  Now 

Q=p-i|PAi  =  i^ili^r'^l  =  i^i~ 


/0(i)  0(1)  0(1) 
0(1)  0(i)  0(1) 

V  0(1)  0(1)  0(1) 


i.e.  Q  is  the  nonpreconditioned  Roe  matrix  which  does  not  have  the  desired  property.  We  therefore  conclude 
that  for  an  upwind  difference  scheme  the  Riemann  solver  should  be  based  on  the  preconditioned  system  and 
not  the  original  scheme.  In  [3]  plots  are  shown  to  illustrate  the  greatly  improved  accuracy  for  low  Mach 
number  flows  when  the  Riemann  solver  is  based  on  the  preconditioning.  Characteristics  in  the  boundary 
conditions  these  should  be  based  on  the  characteristics  of  the  modified  system  and  not  the  physical  system. 
Preconditioning  is  even  more  important  when  using  multigrid  than  with  an  explicit  scheme.  With  the  original 
system  the  disparity  of  the  eigenvalues  greatly  affects  the  smoothing  rates  of  the  slow  components  and  so 
slows  down  the  multigrid  method,  [6]. 

We  conclude  from  the  above  remarks  that  the  steady  state  solution  of  the  preconditioned  system  may 
be  different  from  that  of  the  physical  system.  Thus,  on  the  finite  difference  level  the  preconditioning  can 
improve  the  accuracy  as  well  as  the  convergence  rate. 

Algorithm 


In  terms  of  the  primitive  variables  the  preconditioning  we  consider  is: 


( 

a" 

6^ 

0 

0 

0 

\ 

ocau 

~W 

1 

0 

0 

aav 

(3^ 

0 

1 

0 

0 

0 

0 

1 

(  u 

a 

0 

0 

\ 

a 

u 

0 

0 

+ 

0 

0 

u 

0 

0 

0 

0 

u 

/ 

i  V 

0 

a 

0 

\ 

0 

V 

0 

0 

a 

0 

V 

0 

VO 

0 

0 

V 

/ 

/  ^  \ 

du 

dv 

\dS  J, 

/  \ 

pa 

du 

dv 

\dS 

i  \ 

du 

dv 

XdsJ, 


The  nonpreconditioned  case  corresponds  to  (3^  =  a  =  0.  Let  q  =  ucoi  +  vuj2j  then  the  eigenvalues  of  PD 
are  given  by 

do  =  q  (double) 


=  I  [(1  -a  +  j3^/a^)± 


(8) 


(1  —  a  +  l3'^/a?Y  +  '^(‘^1  +  ‘^2 
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For  general  curvilinear  coordinates,  in  the  ”i”  direction  wi  =  t/^^,  a;2  =  — The  time  step  is  bounded  by 
VOL 

Our  ultimate  goal  is  to  have  a  compressible  code  that  solves  the  incompressible  equations  when  the  input 
Mach  number  is  zero.  So  we  wish  to  use  variables  that  give  us  the  same  result  as  an  incompressible  code  on 
all  levels  of  the  algorithm,  e.g.  flux  computation,  boundary  conditions,acceleration  techniques,  etc.  Hence, 
we  choose  as  our  basic  variables 


W,= 


F  = 


(  P'  \ 

pu 
pv 

V  E'  ) 

(  pu  \ 
pu"^  +  p' 
puv 
\  pH'u  ) 


Q  = 


(  P'  \ 

pu 
pv 

\H'  J 


G  = 


pv 
puv 

pv^  +  p' 

V  pH'v  ) 


dW, 


p  _ 


dt 

where 


=  Pi 


dx  dy  \ 


P  =  P- Poo 


hci 


7-1 


p(u^  +  v"^ 


E'  =  Cpp(T-T’co)-(p-Poo)  + 

=  E->r  Poo -  hoop 
pH'  =  E'  +p'  =  E  +  p-hoop 

We  subtract  the  constants  to  keep  the  quantities  in  scale.  Density  is  now  calculated  from  the  pressure  and 
total  energy.  Because  the  modified  energy  E'  also  contains  the  density  we  get  a  quadratic  equation  for 
the  density.  Choosing  the  positive  square  root  guarantees  that  the  density  is  always  positive.  The  residual 
smoothing  and  multigrid  are  applied  to  p'  and  E'  rather  than  p  and  E.  Thus,  we  duplicate  the  treatment  of 
the  variables  in  a  pseudo-compressible  incompressible  code. 


Pp  =  H-A- 


/  1-i 

—  uB2 
—vB2 

V  ~E4  c+s; 


G+Aoo 

u^B-> 

G-\-hoo 

UvB‘2 

G-\-hca 

uBa 


G^hoo 

UvB‘2 

G+hoo 

V^B‘2 

G+/I00 

vBa 


G^hoo 
—  uB‘2 
G-\-hoo 

G-{-ht3o 


\ 


) 


<?+Aco  G^h 

where  h  =  CpT  =  G  —  A  =  In  the  appendix  we  derive  this  form  of  P 

B.  1  ‘  > 

B2 
Ba 


p- 


/?2  (7  -  l)h  ^2 

T..  a 

+  I2 


/?2 

=  BiH'  + 


a(u^  -I-  v^) 
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We  choose  /?^  =  min  {max  -h  ,a^}  where  prnin  should  have  the  units  of  speed.  The 

choice  of  fimin  is  discussed  in  the  result  section.  In  all  cases  a  =  min  1?  •  We  can  evaluate  this 

efficiently  by  defining  5  =  A  •  (^dp  —  ^  ~  •  Then 

^Pnew  ^ 

^(^P'^^new  —  ^(^P'^'jorig  S^uS 

d{pv)  new  —  di^pV^Q'rig  B2VS 

^^new  ~  dEorig  —  B4S 

These  equations  are  given  for  the  nondimensionalized  variables.  The  nondimensionalization  affects  the 
convergence.  In  some  codes,  p  and  p  are  fixed  in  the  far  field.  This  implies  that  the  speed  of  sound,  a,  is 
also  bounded.  As  the  Mach  number  goes  to  zero  the  pressure  remains  of  order  1  while  the  velocities  go  to 
zero.  Alternatively,  one  can  nondimensionalize  so  that  the  velocities  are  of  order  1  in  the  far  field  and  then 
the  pressure  and  speed  of  sound  go  to  infinity,  unless  one  subtracts  an  appropriate  constant, 

The  boundary  conditions  at  the  far  field  boundary,  for  subsonic  flow,  are  based  on  the  one  dimensional 
theory  of  characteristics  in  the  direction  normal  to  the  boundary.  The  preconditioning  changes  the  form  of 
these  characteristic  variables.  In  differential  form  they  are  given  by 

1  ^  0^ 

_y(„(l_a  +  ^))2+4(i_^)^2^  dp' 

+yKl-a  +  ^))2+4(l-:^)/?2j  dp' 

where  u  is  the  component  of  the  velocity  normal  to  the  boundary.  If  we  consider  low  Mach  numbers  then 
we  can  approximate  these  by 

_  1  dp'  r,  j  dp' 

=  du-\’  — ,  R2  =  du - - 

P0  P0 

which  is  the  same  as  for  the  incompressible  case.  Hence,  at  inflow  i2i,  v  (tangential  velocity)  and  S  are 
specified  while  R2  is  extrapolated  from  the  interior.  We  then  calculate  u  (normal  velocity)  and  the  pressure 
from  Ri  and  R2  and  then  the  density  and  total  energy.  At  outflow  the  role  of  specified  and  extrapolated 
quantities  is  reversed.  At  solid  boundaries  the  normal  momentum  equation  is  used  which  is  not  affected  by 
the  preconditioning. 

Computational  Results 

The  solution  is  advanced  by  a  explicit  Runge^Kutta  method  ([5], [8])  with  residual  smoothing  and  multigrid 
and  no  enthalpy  damping.  In  all  cases  three  levels  of  FMG  multigrid  were  used  with  50  Runge-Kutta  cycles 
on  the  coarser  grids.  Hence,  all  plots  show  the  convergence  for  two  sets  of  50  cycles  and  then  the  convergence 
on  the  finest  mesh.  The  plots  are  of  the  convergence  rate  of  the  residual  of  the  continuity  equation.  For  the 
original  code  this  was  updated  for  the  density  while  in  the  preconditioned  code  it  is  updated  for  the  pressure. 
Nevertheless,  in  the  steady  state  the  residual  of  the  continuity  equation  should  be  the  same  except  for  the 
change  in  the  artificial  viscosity  between  the  two  algorithms.  All  cases  were  run  with  a  matrix  viscosity. 

We  first  present  two  calculations  for  inviscid  flow  about  a  NAG  A  0012.  We  use  a  224  x  32  C  mesh 
and  three  levels  of  multigrid.  The  first  calculation  is  for  inflow  conditions  M  =  0.01,  a  =  1.25^  .  In  this 
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case  we  see  that  the  residual  asymptotes  without  the  use  of  preconditioning  and  that  the  preconditioning 
dramatically  increases  the  rate  of  convergence.  The  use  of  the  preconditioning  adds  only  a  few  percent 
to  the  total  computational  time.  In  the  second  case  we  consider  the  same  geometry  but  with  an  inflow 
of  M  =  0.7,  a  =  1.25^  We  have  also  done  M  =  0.8,  a  =  1.25^  which  results  in  a  minor  slowing  of  the 
convergence  rate.  The  preconditioned  residual  is  the  dotted  line  and  the  original  code  is  the  solid  line. 
Different  parameters  for  the  time  step  and  residual  smoothing  are  needed  with  and  without  preconditioning. 
For  inviscid  cases  we  can  choose  (Smin  as  zero  while  for  the  viscous  cases  ^rnin  =  For  the 

transonic  cases  the  lift  and  drag  coefficients  are  changed  only  minimally  by  the  preconditioning. 

We  next  consider  viscous  flow  about  a  RAE2822  airfoil  on  a  320  x  64  C  mesh  and  5  levels  of  multigrid  on 
the  finest  grid  with  Moo  =  0-01,  a  =  2,79®  using  a  Baldwin-Lomax  turbulence  model  with  Re  =  6,5  million. 
The  residual  history  is  presented  in  figure  3.  Again  the  standard  code  converges  very  slowly  for  these  low 
Mach  numbers.  In  figure  4  we  present  both  the  preconditioned  residual  (dashed  line)  and  the  original  code 
(solid  line)  for  the  same  case  but  Moo  =  0.73.  For  viscous  cases  we  choose  pmin  =  0.4.  Again,  for  the 
transonic  cases  the  lift  and  drag  are  changed  by  about  2  percent  by  the  preconditioning.  For  the  very  low 
Mach  numbers  the  lift  and  drag  coefficients  never  converged  for  the  non-preconditioned  algorithm  and  seem 
to  have  significant  errors.  The  preconditioned  code  gives  much  better  agrees  for  lift  and  drag  for  low  Mach 
numbers. 

We  conclude  with  a  three  dimensional  catse,  inviscid  flow  about  an  ONERA  wing.  In  figure  5  we  display 
the  convergence  rate  for  the  continuity  equation  (normalized  by  the  initial  residual)  for  Mach  numbers  .10, 
.05  and  .01.  We  see  that  the  convergence  rate  is  independent  of  the  inflow  Mach  number.  In  figure  6  we 
plot  the  lift  coefficient  for  the  same  case.  We  again  see  that  the  lift  coefficient  is  essentially  independent 
of  the  Mach  number  except  for  some  slight  compressibility  effects.  Without  preconditioning  there  are  large 
variations  in  the  lift  for  this  set  of  Mach  numbers. 
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Appendix 


To  find  Pp  we  begin  with  the  preconditioner  P5  for  the  variables  dWs  =  {^,du,dv,dSy ^  with  dS  — 
dp  —  a^dp.  We  then  transform  to  dWs  =  (dp,  du^  dv,  dSf  by  multiplying  all  elements  in  the  first  row  of  the 
matrix  by  pa  and  every  element  in  the  first  column  by  This  gives 
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We  then  transform  to  the  conservation  variables  Wc 

r  =  7-  1 


{p,pu,pv,Ey.  This  is  given  by  dWc  =  TidWs-  Let 
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where  G  =  ^  Y  and  a  is  the  speed  of  sound.  This  gives  the  preconditioner  in  conservation  variables. 

Let  Q.  =  ,  Q,  ^  ^  £  _  i  +  -  n).  Then  in 


conservation  variables  P^  =  TiP 

STI^ 

l  +  GQi 

-uQi 

-vQi 

Qi  \ 

UGQ2 

1  U^Q2 

—UVQ2 

uQ2 

Mr  c  — 

vGQ2 

—UVQ2 

1  - 

vQ2 

GR 

—uR 

—vR 

1  +  R  / 

We  next  change  from  Wc  =  {p,  pu^  p,  E)  variables  to  u;'  =  (p, /?u, /?,  E")  variables,  =  E  —  phoo  -hPooj 
dW'  =  T2dWc. 


T2- 


P'  — 


/  1 
0 
0 

\  -ho 


0  0  1  \ 
1  0  0 
0  1  0 
0  0  1/ 


/  1  +  GQi  —uQi 

uGQ2  1  -  U^Q2 
VGQ2  -UVQ2 
\  RG  -uR 


-vQi 
—UVQ2 
1  -  V^Q2 
—vR 


-Qi  \ 
uQ2 
vQ2 
1  +  R  / 


Then  P'c  =  T2(TiPsT£l)T2^ 
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We  finally  change  only  the  time  derivatives  to  the  variables  Wp  by  dWp  —  TzdW^. 


(  r((j  +  /loo) 

T-  ° 

J3-  Q 

\  0 


T3-1  = 


r(G+ftoo) 

0  1 

0  0 

0  0 


To  summarize,  we  begin  with 


-idWs 


oy 


and  transform  to 

We  then  transform  ,  in  conservation  form,  to  the  prime  variables  where  E  is  replaced  by  E'  =  E  —  h^op. 
Finally  we  then  have  that  =  P'  =  T2T1PS  ^  “  ^3^27iP5^i  ^"2  ' 


Finally  we  then  have  that  P^ 
Thus, 


orPp^TaP'e 


^uAB2  1+^ 
^vAB2 
-AJ54 


§1 

K 

^u(3^B2 

K 

-vP^B2 

1-^ 


where  all  quantities  were  defined  in  the  text. 

We  next  show  how  to  convert  any  preconditioner  given  in  streamline  coordinates  and  {■^,du,dv,dS) 
coordinates  to  conservative  variables  in  Cartesian  (not  streamwise)  coordinates.  We  shall  do  this  in  two 
dimensions  but  the  extension  to  three  dimensions  is  straightforward.  Assume  we  are  given  a  preconditioner 
in  streamline  coordinates  and  (^,du,dv,dS)  coordinates  Ps  given  by 

/  Pii  P12  0  0  \ 

p  _  P21  P22  0  0 

-  0  0  P33  0 

VO  0  P55  / 

We  define  rotation  matrices  U,  U~^  to  get  P  in  Cartesian  coordinates. 


10  0  0 
0  cos9  sinO  0 
0  —sinO  cosO  0 
0  0  0  1 


10  0  0 
0  cos9  -’Sin9  0 
0  sin9  cos9  0 
0  0  0  1 


Let  To  get  the  streamwise  direction  we  shall  choose 


cos9  = 


sin9  = 


-h 
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Then  the  preconditioner  in  Cartesian  coordinates  is  given  by  Pear  =  U  and 

/  Pll  Pl2f  Pl2f  0  \ 

,  _  P2lf  P22$  +  PSZ$  (P22-P33)f  0 

P2lf  (P22-P33)fl  P22f^  +  P33fj  0 

Vo  0  0  P55  / 


^  Qii 

Qi2 

<5 13 

0  \ 

Q21 

Q22 

Q23 

0 

Q31 

Qz2 

Qzz 

0 

0 

0 

0 

Qzz  ) 

We  next  introduce  conservative  variables  Wc  as  given  in  the  appendix  by  the  transformation  Ti-  The 
preconditioner  for  conservative  variables  is  then  given  by  Pc  =  Ti  Pear'll  ^ 

We  now  define  the  following  quantities 


Yi 

uQzi 

+  vQzi  +  U'Q41 

Yz 

= 

UQ22  +  vQs2  +  '^Qa2 

Yz 

uQ2Z  +  vQzz  d-  wQ^z 

Ya 

uQ2A  “b  ^^34  "b 

Yz 

•  = 

(7- 

1)(Q55  -  <3ll) 

G 

- 

(t- 

V 

l)g2 

> 

L 

= 

(t- 

!)/> 

Zu 

= 

1 

pa  2 

WYz  Til 

2  p\ 

+  Qss 

Z\2 

= 

1  ^  ,  Qi2 

— 2  uYs  + - 

pa?  p 

Z\z 

1  ^  ,  Qiz 

pa^  p 

Z\A 

^ 

pa2  p 

Ziz 

= 

-^Yz 

pa^ 

Zz^ 

HZn  +  GYi  -  uY2  -  vYz  -  wYa 

Zz2 

= 

HZ12  —  (7  ~  l)pwTi  +Y2  —  uQzz 

Zzz 

HZiz  -  (7  -  Y)P'vYi  +Yz-  uQzz 

Zz4 

= 

HZia  -  (7  -  l)pwYi  +Y4-  uQzz 

Zzz 

HZiz  +  (7  -  PjpQzz 
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Then 


Pc  = 


/  2^11 

uZii  +  GQ2X  “  ^2 
vZii  +  GQzi  -  V3 
+  GQ41  —  y4 
\  ^51 


^12 

^^12  “  i?tzQ21  +  Q22 
^■^12  “  -RwQsi  +  Q32 
'wZi2  —  RuQ^^i  +  Q42 
^52 


^13 

ti^l3  —  RvQ21  +  Q23 
vZia  —  RvQsi  +  Q33 
It; 2^23  ■”  RvQ^i  +  Q43 

Zs3 


^14 

uZi4.  —  R'wQ2i  +  Q24 

vZi4  —  HtoQsi  +  Q34 
iy^i4  —  RwQ^i  +  Q44 
-^54 


^15  \ 

uZis  +  RQ21 

vZis  +  HQ31 

^55  / 
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50  100  150  200  250 

Figure  1:  Convergence  rate  for  inviscid  flow  about  a  NACA0012  with 


NACA0012  M=0.7  a=1^5 


Figure  2:  Same  as  figure  1  with  Moo  =  0 


.70  am 


Figure  3:  Convergence  rate  for  viscous  flow  about  a  RAE2822  airfoil  with  Moo  —  0.01  and  a  —  2.79®  Solid 
line  is  original  algorithm  and  dashed  line  is  the  preconditioned  scheme 


Figure  4:  Convergence  rate  for  viscous  flow  about  a  RAE2822  airfoil  with  Moo  =  0.73  and  a  =  2.79®.  Solid 
line  is  original  algorithm  and  dashed  line  is  the  preconditioned  scheme 
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LHl  Coefficient 


Figure  5:  Convergence  rate  for  inviscid  flow  about  ONERA  wing.  Moo  =  .10,  .05,  .01,  a  =  3.06^ 
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Figure  6:  Lift  coefficient  for  inviscid  flow  about  a  ONERA  wing. 


13 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


1.  AGENCY  USE  ONLYfleave6/an*;|2.  REPORT  DATE  1 3.  REPORT  TYPE  AND  DATES  COVERED 

I  May  1995 _  |  Contractor  Report 

4.  TITLE  AND  SUBTITLE -  FUNDING  NUMBERS 

PRESSURE  UPDATING  METHODS  FOR  THE  STEADY-STATE 

FLUID  EQUATIONS  WU  505-90-52-01 

6.  AUTHOR(S) 

A.  Fiterman 
E.  Turkel 

V.  Vatsa  _  — ^ 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Institute  for  Computer  Applications  in  Science 
and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 
Hampton,  VA  23681-0001 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

National  Aeronautics  and  Space  Administration 
Langley  Research  Center 
Hampton,  VA  23681-0001 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ICASE  Report  No.  95-40 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

NASA  CR-198163 
ICA.SE  Report  No.  95-40 


11.  SUPPLEMENTARY  NOTES 

Langley  Technical  Monitor:  Dennis  M.  Bushnell 
Final  Report 

To  appear  in  the  Proc.  of  the  AIAA  CFD  Conference 

12a.  DISTRIBUTION/ AVAILABILITY  STATEMENT 

U  nclassihed-U  nlimited 


1 12b.  DISTRIBUTION  CODE 


Subject  Category  64 


13.  ABSTRACT  ^  i  ^  i  x  -r 

We  consider  the  steady  state  equations  for  a  compressible  fluid.  Since  we  wish  to  solve  for  a  range  of  spe^s  we 

must  consider  the  equations  in  conservation  form.  For  transonic  speeds  these  equations  are  of  mixed  type.  Hence, 
the  usual  approach  is  to  add  time  derivatives  to  the  steady  state  equations  and  then  march  these  equations  m  time. 
One  then  adds  a  time  derivative  of  the  density  to  the  continuity  equation,  a  derivative  of  the  momentum  to  the 
momentum  equation  and  a  derivative  of  the  total  energy  to  the  energy  equation.  This  choice  is  dictated  by  the 
time  consistent  equations.  However,  since  we  are  only  interested  in  the  steady  state  this  is  not  necessary.  Ihus 
we  shall  consider  the  possibilty  of  adding  a  time  derivative  of  the  pressure  to  the  continuity  equation  and  similar 
modifications  for  the  energy  equation.  This  can  then  be  generalized  to  adding  combinations  of  tune  derivatives  to 
each  equation  since  these  vanish  in  the  steady  state.  When  using  acceleration  techniques  such  ^  residual  smoothing, 
multigrid,  etc.  these  are  applied  to  the  pressure  rather  than  the  density.  Hence,  the  code  duplicates  the  behavior  of 
the  incompressible  equations  for  low  speeds. 


14.  SUBJECT  TERMS 

Preconditioning;  Low  Mach  Flows;  Euler  Equations 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

Unclassified 


NSN  7540-01-280-5500 


15.  NUMBER  OF  PAGES 

15 


16.  PRICE  CODE 

A03  _ 


18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION 
OF  THIS  PAGE  OF  ABSTRACT  OF  ABSTRACT 

U  nclassified  _ _  _ 

Prescribed  by  ANSI  Std.  Z39-18 


